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A Numerical Study of Resistivity and Hall Effects for 
a Compressible MHD Model 

B. Sjogreen * H. C. Yee * 


The effect of resistive, Hall, and viscous terms on the flow structure compared with 
compressible ideal MHD is studied numerically for a one-fluid non-ideal MHD model. The 
goal of the present study is to shed some light on the emerging area of non-ideal MHD 
modeling and simulation. Numerical experiments are performed on a hypersonic blunt 
body flow with future application to plasma aerodynamics flow control in reentry vehicles. 
Numerical experiments are also performed on a magnetized time-developing mixing layer 
with possible application to magnetic/turbulence mixing. 


I. Motivation 

The role various components of the generalized Ohm’s law can play, in many of the physical phenomena 
encountered in physics and engineering is not fully understood. Four of the emerging areas of non-ideal 
MHD modeling are: (a) Plasma aerodynamic flow control for possible drag reduction, normal force control 
and heat transfer minimization of highly maneuverable high speed combat aircraft and reentry vehicles, 
(b) Reentry vehicles at hypersonic speed that develop ionized plasma that creates a magnetic field which 
can affect the design of heat shield, flight control and communication apparatus, (c) Plasma detachment 
or magnetic reconnection (MR) in the design of Tokamak reactors, in plasma rocket exhaust and in the 
dynamics of solar coronal mass ejections and their effect on the structure of the solar wind, the heliospheric 
magnetic field and planetary magnetospheres, and (d) Laser, fusion and high energy density physics. 

Pertaining to areas (a) and (b), for example, in hypersonic blunt body gas dynamics flow, it is known 
that different equations of state and evolutionary equation sets (e.g., non-equilibrium flows) can affect the 
shock standoff distance and bow shock shape, and consequently heating of the blunt body. Plasma injection 
in hypersonic reentry vehicles for possible heat transfer minimization is important in the design of heat shield 
of the vehicle. Pertaining to areas (c) and (d), aside from dealing with general equations of state, multi- 
fluid and multi-phase flows, for some MHD modeling, in order for plasma detachment or MR to proceed, 
there must be some mechanism which can violate the ideal MHD condition in the diffusion region and allow 
plasma to diffuse across the magnetic field. One of the candidates for breaking the ideal MHD condition 
is resistivity. 6, 11 Another candidate is the well known Hall effect and thermal inertia terms. 7,8 A third 
candidate is through turbulent reconnection and reconnection in chaotic magnetic fields. 2 ’® 

In' order to shed some fight on the effect of the resistive and Hall coefficients on the flow structure, 
we first perform numerical study on a one-fluid non-ideal MHD model. Numerical experiments on (a) a 
Mach 15 and 1.5 blunt body flow and (b) a magnetized version of a vortex pairing in the time-developing 
mixing layer problem that was previously investigated by the authors and collaborators 16 are conducted. 
For the blunt body problem the effect of the resistive and Hall coefficients on the shock standoff distance 
and bow shock shape compared with the perfect gas and ideal MHD cases will be illustrated. The effect 
of the structure of the bow shock as a function of the plasma j3 (/? p ) will also be investigated. Due to the 
simple flow structure of the blunt body problem, only second-order shock-capturing methods are used for 
the simulation. The magnetized version of a vortex pairing in the time-developing mixing layer problem is 
a time-accurate computation with complex flow structure. Our adaptive numerical dissipation control in 
high order filter schemes 12,13,16, 17 > 19 ~ 21 i s used for the simulation. It was shown that our scheme is stable 
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and accurate for a wide spectrum of flow physics ranging from long time wave propagation of smooth flows 
to high speed multiscale turbulent flows including strong shock waves. Two unique features of our scheme 
are (a) Solving the ideal conservative MHD system (a non-strictly hyperbolic system of conservation laws) 
without having to deal with its incomplete eigensystem, and (b) Adaptive numerical dissipation control that 
permits, at ease, the minimization of the divergence of the magnetic field (V-B) numerical error. Since 
second-order shock-capturing schemes are well known and can be found in many text books and review 
articles, 15 the next section gives only a brief review of our high order filter scheme. Numerical experiments 
on the one-fluid resistive MHD model will be given in Section III. 

II. Designing Adaptive Low Dissipative High Order Schemes for the 

Compressible MHD Equations 

Our adaptive numerical dissipation control in high order filter schemes 17 consists of automatic detection 
of different flow features by distinct sensors to signal the appropriate type and amount of numerical dissi- 
pation/filter where needed while leaving the rest of the region free of numerical dissipation contamination. 
These scheme-independent detectors are capable of distinguishing shocks/shears, flame sheets, turbulent 
fluctuations and spurious high-frequency oscillations. In addition, these sensors are readily available as an 
improvement over existing grid adaptation indicators. The detection algorithm is based on an artificial com- 
pression method (ACM) (for shocks/shears), and redundant multiresolution wavelets (WAV) (for all types 
of flow feature). The ACM and wavelet filter schemes with sixth-order spatial central scheme for both the 
inviscid and viscous flux derivatives and a fourth-order Runge-Kutta method (RK4) are denoted by ACM66 
and WAV66. A second-order MUSCL spatial scheme with a second-order Runge-Kutta method (RK2) is 
denoted by (MUSCL). 

In our previous work, extensive grid convergence studies were performed using WAV66 and ACM66 for 
many well tested model problems. In many instances, grid convergence was achieved by ACM66 and WAV66 
but not by MUSCL using the same grid sequence. More accurate solutions were obtained with WAV66 
and ACM66 than with MUSCL, which requires similar CPU time, and a fifth order weighted ENO scheme 
(WEN05-RK4), which is nearly three times as expensive. See 12, 13 > 1 6.i7, 19-21 f or details. The following 
subsections give an brief overview of features (a) and (b) and a description of the filter method. 

A. Solving Conservative MHD Systems Using Symmetrizable Eigenvectors 

Consider the 3-D conservative and symmetrizable 5 (non-conservative) forms of the ideal compressible MHD 
equations, 

U t + V • F = 0 (conservative), (1) 

U t + V • F = S (symmetrizable), (2) 

f P 
pu 
pv 

U = pw 
e 

B x 
By 
\B Z 

The velocity vector u = ( u,v,w) T , the magnetic field vector B = (B x , B y , B Z ) T , p is the density, and e is 
the total energy. The notation B 2 * = B 2 + B 2 + B 2 is used. The pressure is related to the other variables by 

P = (7 ~ l)(e - \p{u 2 + v 2 + w 2 ) ~ \(Bl + B\ + B 2 )). 

For plasmas and monatomic gases, 7 = 5/3. The vector on the right hand side of ( 2 ) is the non-conservative 
porti on of t he sym m etri zabl e MHD e quati ons an d is frequently referred to in the literature as a source term 
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vector. The conservative and symmetrizable forms of the non-ideal compressible MHD 3 (viscous, resistive 
and Hall MHD) take the form 

17 t + V ■ F = F*,, 

+ V • F = F„ + S', 

1 T 
F v = [ 0 divr f v5 ~(AB - VdivB) - /? h V x ((V x B) x B) ]. . 

<j 

The fifth component of F„ is 

/„ 5 = div(u T r) + divh — — div((V x B) x B) - /3/idiv (((V x B) x B) x B) . 

u 

The vector F„ includes viscosity, resistivity, and Hall effect with r being the viscous stress tensor, a the 
conductivity coefficient, flh the strength of the Hall effect (fi with subscript “h”), and h the heat flux. The 
notation for the plasma j3 (/ 3 with subscript “p”) /3 P = plasma pressure / magnetic pressure will be used. 

An important ingredient in our high order filter method is the use of the dissipative portion of high- 
resolution shock-capturing schemes as part of the nonlinear filters for discontinuity capturing. These non- 
linear filters usually involve the use of approximate Riemann solvers (knowledge of the eigenstructure of the 
system). Due to the fact that the conservative MHD system (1) is a non-strictly hyperbolic conservation 
law, we proposed in 13,20 to solve (1) using the eigenvectors of the symmetrizable system (2). See 13,20 for 
details. 

B. Divergence- Free Preserving High Order Filter Methods 

The filter method consists.of two steps, a divergence-free preserving (base scheme) step (not involving the use 
of approximate Riemann solvers or flux limiters) and a filter step (usually involving the use of approximate 
Riemann solvers and flux limiters). The type of base schemes is very general. The adaptive filter is scheme 
independent and can be used in conjunction with spectral, compact or non-compact spatially central base 
schemes. In order to have good sho ck-capturing capability, improved nonlinear stability related to spurious 
high frequency oscillations, and no interference with the accuracy of turbulent fluctuation, the blending of a 
high order nonlinear filter and a high order linear filter was proposed by Yee & Sjogreen . 17,20 The nonlinear 
filter consists of the product of an ACM or wavelet flow sensor and the nonlinear dissipative portion of 
a high- resolution shock-capturing scheme. The high order linear filter consists of the product of another 
sensor, a tuning parameter and a high order centered linear dissipative operator that is compatible with the 
order of the base scheme being used. Three forms of the nonlinear dissipative portion of shock-capturing are 
considered, namely: 

• Dissipative portion of the fifth-order WENO scheme (WEN05). It can be obtained, e.g., in the in- 
direction by taking the full WEN05 scheme in the ^-direction and subtracting D^Fj (sixth-order 
central differencing). 

• Dissipative portion of the a second-order MUSCL scheme . 16 

• Dissipative portion of the Harten-Yee TVD scheme . 16,20 

The resolution of the above forms of the dissipative portion of shock-capturing schemes are comparable and 
work well on all of our numerical tests for both the ideal and resistive viscous MHD systems. See Yee & 
Sjogreen 21 and references cited therein for details. 

The nonlinear filter we proposed, if applied to the entire MHD system, will not preserve the divergence 
free magnetic field condition in general. For the computations in this paper, the “No filter on B” option 
is chosen. That is, the nonlinear filter step only applies to the first five equations of (1) or (2). With the 
divergence free spatialbase scheinegthe divergence free property should be preserved for periodic boundary- 
conditions on uniform Cartesian grids. Extensive grid convergence comparison of the “no filter on B” with 
the “filter all of the MHD equations” options were presented in . 20 
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III. Numerical Experiment 

For the numerical experiment, we chose (a) a Mach 15 blunt body and the De Sterck and Poedts 14 (Mach 
1.5, — 0.4) blunt body flows, and (b) a magnetized version of a vortex pairing in time-developing mixing 

layer problem that was previously investigated by the authors and collaborators. 12,16,17 
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Fig. 1 Blunt body problem with M = 15: Density contours comparison: Perfect Gas, Ideal MHD 0 P = 5, 
Ideal MHD /3 P = 10, Resistive MHD Re = 1000, a = 100, 0 p = 5,0 h = 0, IfosMiw, MHD Re = 1000, u = 
100, 0 P — 5,(3 h — 0.5. Second-order MUSCL using a 121 x 41 grid 



Fig. 2 Blunt body problem with M = 1.5, 0 P = 0.4: Density contours comparison: Second-order Lax- 
Friedrichs, Harten-Yee and MUSCL using a 81 x 81 grid 


Blunt Body Results (M = 15 & 1.5): For this problem, only second-order shock-capturing schemes using 
a second-order Runge-Kutta temporal discretization (RK2) solving the symmetrizable MHD system (2) are 
employed for the study. For time-marching to the steady states, the schemes have problems with convergence 
if we solve the conservative MHD system (1). The reason is not certain and is currently under investigation. 

Figure 1 shows the comparison among perfect gas, ideal MHD with two plasma 0 (0 P = 5, 10), resistive 
MHD with Re = 1000, a = 100, /3 P — 5, and Hail MHD with. the Hall coefficient 0h = 0.5. The computation 
shown uses a second-order MUSCL scheme with the van Albada limiter using a mildly stretched parabolic 
curvilinear 121 x 41 grid. Grid refinement indicates that the 121 x 41 grid is sufficient to obtain an accurate 
solution. The shock standoff distance and the bow shock shape are very different among the perfect gas, 
ideal MHD and resistive MHD. 

For the Hail coefficient 0h < 0.5, the Hall effect is very small compared with the resistive MHD for 
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Fig . 3 Blunt body problem M = 1.5, fi p = 0.4; Density contours of the unsteady evolution on an 161 x 161 
grid using MUSCL at 9 different time instances : T=80, 120, 160, 200, 560, 900, 1300, 1600, 2000 
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Fig. ^ Blunt body problem with M = 1.5 ,fi p = 0.4: Comparison: Density contours of the steady solutions by 
the second-order Lax-Friedrichs using 81 x 81, 161 x 161 and 321 x 321 grids 
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the same Re = 1000, a = 100 and plasma B p in the range (0.4,10). As the Hall coefficient 0h increases, 
the problem becomes very stiff and the scheme experiences convergence problems. Future study using the 
implicit stiffly stable local time-stepping approach is planned to further investigate the large fih issue. 

For the ideal MHD with large plasma P p in the range of (0.15, 10) there is no noticeable change in the 
shock standoff distance and bow shock shape. For /3 P below 0.15, the bow shock standoff distance and bow 
shock shape change (figures not shown). 

In the supersonic and low hypersonic range, the bow shock shape is more sensitive to the same range of 
P p than the higher Mach number case (e.g., M = 15). To illustrate this point, we simulate the same problem 
as in De Sterck and Poedts. 14 For M = 1.5 and /? p below 0.7, DeSterck-Poedts concluded that intermediate 
shocks would developed using the second-order Lax-Friedrichs scheme with an 81 x 81 coarse grid. Their 
numerical solutions deviate from the standard bow shock shape near the stagnation region. See 10 for an 
overview of stability of intermediate shocks in conjunction with out of plane perturbations. 

We perform an extensive grid refinement study (41 x 41 to 1281 x 1281 grids) on the De Sterck. and 
Poedts blunt body problem for M = 1.5 and /3 P = 0.4 with inflow speed u = 1.5 in both time-accurate 
and time-marching to the steady state approaches. Figures 2-4 show some of the computations. As the 
grid is refined (larger than 101 x 101), there is no conclusive evidence that these intermediate shocks will 
settle down to a definite steady solution, a limit cycle or a spurious solution (numerical artifact). For an 
81 x 81 grid, the second-order Lax-Friedrichs, Harten-Yee and MUSCL all converge to the same steady state 
(time-accurate or time-marching simulations). As we refine the grid, Harten-Yee and MUSCL no longer 
converge to the same steady state solution, but rather have different and very complex iregular structure 
in the vicinity of the stagnation region (figures not shown). Figure 3 shows the time-accurate computation 
by MUSCL with an 161 x 161 grid (similar unsteady behavior for Harten-Yee). We simulate the problem 
with time-marching mode with the same grid, the numerical solution does not settle down but yields similar 
flow structures for MUSCL and Harten-Yee. This is not case with second-order Lax-Friedrichs. Figure 4 
shows the time-accurate computation by the second-order Lax-Friedrichs using three grids. All the solutions 
converge to the same steady state. The time-marching approach using the second-order Lax-Friedrichs also 
converges to the same steady state with the same grid. For a fine grid such as 641 x 641 or larger, even 
the diffusive Lax-Friedrichs scheme does not settle down to the same asymptote (on the large number of 
iterations that we used). Without further research, our preliminary conclusion is that we cannot conclude if 
the asymptoes are stable. 

Vortex Pairing in Time-Developing Mixing Layer Result: Figure 5 shows the schematic of the flow 
condition of the gas dynamics model studies in. 16 For the magnetized case, in addition to the gas dynamics 
initial conditions (ICs) and boundary conditions (BCs) indicated in Fig. 5, the initial magnetic field is 0.1 in 
the ^-direction and zero in the other two directions. A grid mildly stretched in the ^-direction and uniform 
in the ^-direction is used. In Yee and Sjogreen, 22 a study of the resistive effect on the same one-fluid MHD 
model (1) and (2) without the Hall term was performed on this vortex pairing problem. We first summarize 
the result from that study first. 

Figure 6 shows the comparison of our ACM66-RK4 scheme with MUSCL-RK2. ACM66 required at 
least 50% fewer grid points per direction with similar resolution as MUSCL and WEN05-RK4 (figures not 
shown). Figure 7 shows the comparison of ideal MHD with viscous resistive MHD for Re = 10 3 with 5 
different conductivities a = 10 6 , 10 4 , 10 3 , 100 and 50 at time T = 90, The ACM66 solution converges for the 
viscous resistive model for a = 10 6 , 10 s , 10 3 , 100 and 50 using grid sizes 1601 x 1601, 1601 x 1601, 801 x 801, 
201 x 201 and 201 x 201, respectively. For ideal MHD, the fine scale structures are quite well resolved using 
a 1601 x 1601 grid. However, grid convergence has not been obtained for this case. Comparing Figs. 6 and 
7, one can see that lower accuracy and diffusive shock-capturing schemes such as MUSCL for the ideal MHD 
is similar to that of the solution of added Re = 1000 and a = 10 3 by ACM66. See Fig. 6(c) and Fig. 7(d). 
The resolution of WEN05 and Harten-Yee TVD schemes exhibit a similar diffusive pattern as MUSCL. 

Figure 8 shows the comparison of viscous resistive MHD for Re = 10® with 6 conductivity coefficients 
a = 10 6 , 10 4 ,10 3 , 100,50 and 10 at T - 90: The AGM66-solution converges for a = 10 3 , 100 and 50 using a 
grid of 1601 x 1601, 801 x 801 and 401 x 401, respectively. For the three higher a = 10 6 , 10 5 , 10 3 , the fine 
scale structures are quite well resolved using a 1601 x 1601 grid. However, grid convergence has not been 
obtained for these three a values on the very fine scale structures of the flow. 

The study in 22 gives some insights on th e effect of the resistive and viscous terms compared with ideal 
MHD. The most interesting result is that without adaptive numerical dissipation control, the commonly used 
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Vortex Pairing in a Time-Developing Mixing Layer 

(M c =0.8, Re=1000, T R =300K, Prandtl #=0.72) 


Fig. 5. 


Normalized with vorticity thickness: 

Ui - u 3 

" ~ (du/dy) mam 

T & c: determined by assuming constant stagnation enthalpy 
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Inital shear profile: U = 0.5 tanh (2y) 
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Initial perturbations: 
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Vortex pairing in time-developing mixing layer gas dynamics problem. 
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(a) ACMm?101 x 201 (b) ACM66, fc01 x 801 



(c) MUSCLl *201 x 201 (d) MUSCL * 801 x 801 


Fig. 6. Comparison of the. temperature contours of ACM66 (top row) with MUSCL (bottom row) at T — 90 
for ideal MHT) using 201 x 201 and 801 x 801 grids. 
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Fig. 8. Temperature contours of ACM66 for non-ideal MHD with Re = 10® and conductivities a 
10 6 , 10 s , 10 3 , 100, 50, 10 atT = 90. 
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Fig. 9. Temperature contours of ACM66 for three fih for non-ideal MED with Re = 10 5 and conductivities 
a = 10 3 , 10 5 using a 801 x 801 grid atT = 90. 


shock-capturing schemes such as MUSCL, TVD and WEN05 solving the ideal MHD produce solutions as if 
added physical dissipation were present. 

With the Hall term included, for small Hall coefficient @h < 0.2 there is not much effect on the overall 
flow structure over the pure resistive MHD with the same Re and plasma f3 p . As the Hall coefficient is 
increased beyond 0.2, the problem becomes more stiff and the computation using our filter scheme is stable 
only with grids that are smaller than 801 x 801. Although the more diffusive MUSCL and Harten-Yee are 
stable for a denser grid, the resolution is similar to the coarser grid solution by the filter method. Figure 9 
shows three computations with different Hall coefficient (3h and conductivity coefficient. With the studied 
flow parameters, the flow patterns deviate from the Hall MHD slightly. For larger fa, all considered method 
becomes unstable even with smaller re and a values. 

Concluding Remark: Our preliminary study on the two test cases shows the effect of the resistive and 
Hall coefficients on the flow structures compared with the ideal MHD. The result of the two blunt body 
problems indicate plasma injection can alter the shock standoff distance and heating. The study also sheds 
some light on a simplified model related to solar wind physics. Both test cases indicate that the Hall term 
with large Hall coefficient poses a challenge in numerical modeling and simulation. Additional investigation 
is planned to overcome the numerical instability for large Hall coefficients. 
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